% 假设AIS数据存储在一个结构体数组中，每个元素包含Time, Lon, Lat, Course, Speed
clc
clear
currentDir = pwd; % 获取当前目录路径
parentDir = fileparts(currentDir); % 获取当前目录的上一级目录路径
% grandParentDir = fileparts(parentDir); % 获取当前目录的上两级目录路径
% greatGrandParentDir = fileparts(grandParentDir); % 获取当前目录的上三级目录路径

%指定需要查找的文件夹
subdir = 'DATA';
%指定需要查找的文件名
filename = '413168000(NoRepeated).mat';%413168000.mat'; %matlab中包含的是剔除
%构造需要查找的文件的完整路径
full_path = fullfile(parentDir, subdir, filename);

load( full_path)  %此处使用的数据是原始org数据经过后面的[]处理后得到的数据

%% 数据预处理
data(data.Speed<1,:)=[];%锚泊状态不处理
aisdata=data;

% 提取各列
Time = aisdata.time;
Lon = aisdata.Lon;
Lat = aisdata.Lat;
Course = aisdata.Course;
Speed = aisdata.Speed;

% 定义合理的速度范围
maxSpeed = 15; % 假设最大速度为30节
% 假设Speed是以节为单位的速度，需要转换为米/秒
Speed_mps = Speed * 0.514444; % 1节 = 0.514444米/秒

% 计算时间间隔，这里假设Time是以秒为单位的时间戳
time_diff = diff(Time);

% 计算理论上的行驶距离（米）
theoretical_distance = Speed_mps(1:end-1) .* time_diff;

% 设置阈值，这里假设阈值为100米
threshold = 100;

% 计算实际行驶距离（米）
actual_distance = sqrt(diff(Lon).^2 + diff(Lat).^2) .* 111319.9; % 将经纬度差转换为米，假设地球半径为6371公里

% 识别异常点
anomalies = actual_distance > (theoretical_distance + threshold);

% 标记异常点
Lon(anomalies) = NaN;
Lat(anomalies) = NaN;
LL.ais_deleted(:,1)=Lon(anomalies==0);
LL.ais_deleted(:,2)=Lat(anomalies==0);
%计算剔除插值后的轨迹点偏转角的标准差
LL.degrees_org = calculate_track_degree(table2array(aisdata(:,{'Lon','Lat'})));
LL.degrees_deleted = calculate_track_degree(LL.ais_deleted);

LL.std_org=std(LL.degrees_org);
LL.std_deleted = std(LL.degrees_deleted);

% 插值经度和纬度
Lon_interp = interp1(Time, Lon, Time, 'spline');
Lat_interp = interp1(Time, Lat, Time, 'spline');

% 替换异常点
Lon(isnan(Lon)) = Lon_interp(isnan(Lon));
Lat(isnan(Lat)) = Lat_interp(isnan(Lat));



% 分段三次Hermite插值
Lon_pchip = pchip(Time, Lon, Time);
Lat_pchip = pchip(Time, Lat, Time);

% 替换异常点
Lon(isnan(Lon)) = Lon_pchip(isnan(Lon));
Lat(isnan(Lat)) = Lat_pchip(isnan(Lat));

% 准备数据
inputs = [Lon(1:end-1), Lat(1:end-1), Course(1:end-1), Speed(1:end-1)]'; % 输入层
targets = [Lon(2:end), Lat(2:end)]'; % 目标层

% 创建BP神经网络
net = feedforwardnet(10); % 一个隐藏层，10个神经元

% 训练网络
[net, tr] = train(net, inputs, targets);

% 使用网络进行预测
outputs = net(inputs);

% 提取预测的经度和纬度
Lon_pred = outputs(:, 1);
Lat_pred = outputs(:, 2);

% 替换异常点
Lon(isnan(Lon)) = Lon_pred(isnan(Lon));
Lat(isnan(Lat)) = Lat_pred(isnan(Lat));
LL.ais_spline=[Lon Lat];


%计算偏转角标准差,以及计算轨迹曲率
LL.degrees_org = calculate_track_degree(table2array(aisdata(:,{'Lon','Lat'})));
LL.degrees_splined = calculate_track_degree([Lon,Lat]);


LL.std_spline=std(LL.degrees_splined);

LL.mean_spline=mean(LL.degrees_splined);

figure
subplot(2,1,1)
plot(LL.degrees_org,'ko-');
subplot(2,1,2)
plot(LL.degrees_splined,'ko-');